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1 .  PROGRESS  REVIEW 

The  flow  physics  associated  with  strong  viscous-inviscid  interactions 
is  generally  quite  complex.  Computational  methods  to  evaluate  such 
aerodynamic  phenomena  have  typically  been  associated  with  full,  very  time 
consuming,  Navier-Stokes  (NS)  solvers  or  with  approximate  matched  inviscid- 
interacting  boundary  layer  (IBL)  procedures.  The'purpose  of  the  current 
AFOSR  research  program  is  the  development  and  application  of  a  global 
pressure  relaxation  or  flux-split  pressure  differencing  procedure  for 
assessing  steady  and  transient  flows,  respectively. 

A  reduced  Navier-Stokes  (RNS)  system  that  is  more  computationally 
efficient  than  full  NS  solvers  and  more  accurate  and  less  cumbersome  than 
matched  viscous  and  inviscid  (IBL)  methods  has  been  shown  to  apply  to  a 
significant  class  of  aerodynamic  problems.  This  RNS  system  is  a  composite 
of  the  full  Euler  and  boundary  layer  equations  and  is  discretized  to 
optimize  the  numerical  representation  of  viscous  and  inviscid  regions, 
res pecti vely. 

Two  primitive  variable  RNS  systems  have  been  considered.  The  first  is 
for  the  pressure/velocity  variables  and  represents  a  major  upgrading  of  the 
parabolized  Navier-Stokes  or  PNS  methodology.  This  allows  for  full  elliptic 
or  upstream  interactions,  including  flow  separation,  shock  capturing,  etc. 
The  second  is  a  composite  pseudo-potential/ vortical-velocity  model  that 
represents  an  extension  of  full  potential  methodology  to  viscous  interacting 
flows.  This  formulation  allows  for  the  numerical  advantage  of  potential 
flow  methods  and  those  of  boundary  layer  methods  to  be  included  in  a  single 
and  integrated,  composite  RNS  solver. 

The  pressure/ velocity  formulation  has  been  applied  for  two-dimensional 
steady  and  transient  flows  involving  airfoil  viscous/inviscid  interaction, 
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jet/base  flow  interaction,  transonic  and  supersonic  interactions.  The 
composite  velocity  formulation  has  been  applied  for  subcritical  and 
transonic  viscous  airfoil  and  boattail  aerodynamics.  Both  procedures  are 
currently  being  developed  and  applied  for  three-dimensional  flows  over 
afterbodies  and  wing-body  configurations,  where  axial  and  secondary  flow 
separation  and  strong  pressure  interaction  occurs.  Numerous  publications, 
presentations  and  dissertations  have  resulted  from  this  research  activity. 
Recent  papers  and  presentations  are  detailed  in  the  reference  section.  A 
review  of  progress  during  the  year  is  given  in  the  summaries  and  results 
that  follow.  Detailed  information  can  be  found  in  the  references,  most  of 
which  have  been  submitted  for  AFOSR  review.  A  summary  of  highlights  is 
detailed  at  the  end  of  this  report. 

1 . 1  Pseudo-Potential/RNS  Progress  and  Solutions 

This  methodology  has  been  developed  for  the  viscous  RNS/potential  and 
viscous  RNS/Eul er  systems.  The  effect  of  Euler  versus  potential  outer  flow 
conditions  on  shock  capturing,  entropy/ vorti city  production  and  viscous 
interaction,  in  particular  separation,  has  been  evaluated.  A  major  result 
obtained  during  the  investigation  relates  to  the  appropriate  'conservation' 
differencing  that  (i)  accurately  convects  the  vorticity  generated  behind  the 
embedded  shock  and  (ii)  reproduces  the  correct  vorticity/ entropy  across  the 
shock.  A  combined  conservation/ partial  conservation  form  leads  to  the 
capturing  of  very  sharp  shocks  with  accurate  entropy  jumps  and  without  any 
spurious  generation  of  entropy  in  subsonic  regions  or  near  the  leading 
edges.  Recent  solutions  for  the  pressure  and  entropy  distribution  on  a  NACA 
0012  airfoil  are  shown  in  figures  la  and  1b.  As  seen  from  the  figures,  the 
entropy  jump  is  dominant  only  at  the  shock.  The  overshoot  in  entropy  is 
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quite  physical  and  typically  occurs  within  the  structure  of  a  shock.  In  the 
present  case,  shock  capturing  is  determined  by  the  momentum  shift  principle 
of  Osher-’Enquist .  The  effect  of  outer  potential/Euler  conditions  on  the 
separation  bubble  is  given  in  figure  1c.  The  skin  friction  is  extremely 
sensitive  to  the  outer  flow  formulation.  A  detailed  paper  on  this  aspect  of 
the  composite/RNS  investigations  was  presented  at  a  symposium  at  the 
Polytechnic  University  in  August,  1987.  The  full  paper  will  appear  in  a 
special  issue  of  the  journal,  Computers  and  Fluids  in  August  1988. 

This  procedure  has  now  been  extended  to  three  dimensional  flows.  The 
CV/RNS  equations  in  a  generalized  non- orthogonal  coordinate  system  are 
solved  by  the  consistent  coupled  strongly  implicit  (CSIP)  space  marching 
procedure.  The  viscous  analysis  and  early  results  will  be  described  in 
section  1.5.  Application  of  a  three-dimensional  preconditioned  conjugate 
gradient-CSIP  for  the  computation  of  inviscid  flows  past  the  afterbody 
configuration  (figure  2a)  has  been  completed.  Both  subsonic  and  transonic 
solutions  have  been  obtained.  The  pressure  distribution  along  the 
centerline  on  the  top  side  of  the  body  is  depicted  in  figure  2b  for  a  free 
stream  mach  number  of  0.8.  A  comparison  with  a  partial  solution  obtained 
elsewhere  by  a  panel  method  is  also  shown  for  the  boattail  region.  The  mach 
contour  along  the  top  flap  is  shown  In  figures  2c.  The  flow  is  mildly 
supersonic  on  the  boattail. 

A  Ph.D.  student,  Raymond  Gordnier  has  been  visiting  at  the  AFWAL/WPAFB, 
with  Professor  Rubin  during  his  sabbatical  leave  and  now  on  a  weekly  basis, 
to  continue  the  evaluation  of  the  three-dimensional  viscous  afterbody 
interactions.  This  study  should  be  completed  later  this  year,  see  section 


1 . 2  Unsteady  RNS  Flows 

The  RNS  formulation  has  also  been  considered  for  the  computation  of 
unsteady  flow  past  airfoils  at  incidence.  For  this  purpose,  time  consistent 
CSIP,  time  consistent  forms  of  global  and  alternating  direction  line 
relaxation  and  a  sparse  matrix  direct  solver  (DS)  have  been  considered.  It 
should  be  emphasized  that  the  present  RNS  formulation  predicts  the  correct 
viscous-inviscid  interaction  at  the  trailing  edge  and  does  not  require  any 
second-order  upwinding  or  fourth-order  damping  typically  found  in  time 
dependent  Navier-Stokes  solvers.  The  RNS  formulation  also  allows  for  larger 
time  increments  (At)  and  for  more  efficient  placement  of  far  field 
boundaries  and  application  of  boundary  conditions.  A  paper  containing 
additional  details  has  appeared  in  the  AIAA  journal  in  July  1987. 

A  major  result  from  the  unsteady  computations  relates  to  a  hystersis 
phenomena  that  has  been  observed  as  the  airfoil  oscillates  between  0°  and 
13°  incidence.  The  separation,  lift  and  force  patterns  are  quite  different 
for  the  upward  and  downward  motions.  For  example  converged  steady  solutions 
presented  in  figure  3a  are  obtained  at  5.7°  incidence  when  this  angle  is 
approached  from  below.  An  unsteady  solution  is  generated  when  this  angle  is 
approached  with  initial  conditions  generated  from  a  vortical  solution  at  a 
higher  incidence.  The  skin  friction  for  both  the  steady  and  unsteady 
conditions  at  the  same  incidence  is  depicted  in  figure  3a*  Such  phenomena 
have  been  observed  experimentally  but  has  not  been  reported  in  any  earlier 
computations.  Further  analysis  and  grid  resolution  studies  have  validated 
these  results.  Unsteady  flow  solutions  of  the  complete  (ijnm)  form  of  the 
Navier-Stokes  equations  using  the  DS  have  also  been  obtained  for  comparison 
purposes.  These  were  presented  in  the  last  progress  report.  Additional 
information  is  given  in  section  1.6. 


The  lift  coefficient  and  streamline  patterns  at  a  somewhat  larger 
incidence  and  Mach  number  are  shown  in  figures  3b, c.  These  results  are  in 
excellent  agreement  with  full  time-dependent  Navier-Stokes  solutions,  that 
require  considerably  greater  computational  resources  and  times. 

Additional  results  focus  on  an  alternating  direction  predictor- 
corrector  method  for  shedding  phenomena  on  airfoils  and  the  DS  for  transonic 
flows  and  for  laminar  flow  breakdown  on  airfoils  at  incidence.  More 
detailed  solutions  are  contained  in  the  Ph.D.  dissertations  of  S. 
Ramakrishnan  and  E.  Bender  and  will  also  be  presented  in  papers  at  the 
National  Fluid  Dynamics  Conference  this  June  in  Cincinnati. 

1 . 3  Multigrid  Acceleration  of  RNS  Formulations 

In  the  global  pressure  or  pseudo-potential  relaxation  procedures, 
downstream  information  typically  propagates  upstream,  one  mesh  point  during 
each  sweep.  For  flows  with  significant  upstream  influence  and  on  finer 
meshes,  the  rate  of  convergence  will  necessarily  be  significantly  reduced. 

In  order  to  enhance  the  convergence  rate,  alternating  direction  (predictor*' 
corrector)  and  uni-directional  multigrid  procedures  have  been  developed  and 
applied  successfully.  The  latter  is  discussed  here.  Since  the  pressure  or 
"potential”  is  the  only  variable  that  is  being  'relaxed',  a  one-dimensional 
multigrid  procedure,  applied  to  the  PV/RNS  equations,  has  proven  very 
effective  for  highly  stretched  grids  and  separated  flow  regions.  Standard 


multigrid  procedures  fail  in  these  circumstances.  Both  inviscid  and  viscous 
solutions  have  been  computed  to  test  the  applicability  of  this  procedure. 
Transonic  flow  past  a  biconvex  airfoil,  separated  flow  over  a  trough 
configuration  and  trailing  edge  flows  with  large  pressure  interaction  have 
been  evaluated,  as  these  are  representati ve  of  the  interactive  flow 
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phenomena  that  will  arise  in  the  more  complex  configurations  under 
consideration.  In  each  case,  the  multigrid  procedure  accurately  captures 
the  physical  phenomena  at  a  significantly  reduced  computational  cost, 
figure  4a.  In  most  cases  the  standard  multi  grid  approach  either  failed  or 
provided  marginal  improvements  in  the  convergence  rate.  For  a  deep  trough 
(never  considered  by  previous  investigators),  separated  flow  solutions  were 
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obtained  for  a  moderate  Reynolds  number  (8x10  ).  A  fine  grid  solution  for  a 

moderate  trough  is  depicted  in  figure  4a.  On  the  deep  trough,  for  Re  =  10^, 
for  a  very  fine  mesh,  a  breakdown  of  the  laminar  flow  solution  in  the 
separated  flow  region  is  observed.  As  seen  from  figure  4c,  the 
recirculation  bubble  becomes  locally  unstable  and  the  computation  diverges 
locally  (not  globally).  A  similar  phenomena  has  also  been  observed 
previously  by  the  present  investigators  for  the  sine-wave  airfoil  (see 
section  1.6).  This  is  indicative  of  an  apparent  laminar  flow  instability 
that  is  seen  numerically  on  the  fine  meshes  possible  with  the  multigrid 
procedure,  e.g.  500  mesh  points  across  the  trough.  This  instability  may 
signal  inadequate  resolution,  e.g.,  a  small  secondary  vortex  forms,  or 
incipient  laminar  flow  transition,  or  the  occurrence  of  dynamic  stall,  which 
has  also  been  predicted  by  large  Re  asymptotic  methods.  The  present 
results,  which  were  presented  at  the  AIAA  meeting  in  Honolulu,  June  1987, 
provide  a  severe  warning  for  coarse  mesh  calculations.  Recent  results 
obtained  with  full  NS  solvers  tend  to  confirm  the  RNS  predictions  and  lend 
credence  to  the  RNS  approximation  even  for  these  flow  conditions. 

1 .4  Solution  of  RNS  Equations  at  High  Supersonic  Mach  Numbers 

The  global  pressure  relaxation  has  also  been  used  for  the  solution  of 
high  supersonic  flows  past  sharp  nosed  bodies.  For  mach  number?  greater 
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than  3»  all  bow  shock  waves  are  fitted  and  imbedded  shocks  are  captured. 


Two  axisymmetric  configurations  have  been  considered,  i.e.,  a  cone-cylinder 


boattail  and  a  simple  forebody  configuration  of  a  fighter  aircraft.  The 


calculated  mach  contours,  pressure  coefficient  and  skin  friction  for  the 


forebody  are  shown  in  figures  5a, b.  At  the  higher  mach  numbers,  the 


separated  flow  region  that  forms  on  the  boattail,  due  to  shock- boundary 


layer  interaction,  weakens  and  the  shock  becomes  highly  oblique.  The  shock 


that  develops  near  the  canopy  interacts  directly  with  the  bow  shock.  This 


shock-shock  interaction  is  captured  reasonably  by  the  present  technique. 


The  present  formulation  is  highly  efficient  for  large  mach  number  flows. 


This  is  evident  from  figure  5c  where  the  rates  of  convergence  for  the  global 


relaxation  procedure  are  presented.  The  method  converges  very  quickly,  even 


without  the  multigrid  accelerator.  Additional  results  for  an  axisymmetric 


surface  cavity  have  also  been  obtained,  see  figure  5d ,  and  turbulent  flow 


computations  have  been  initiated. 


1 .5  Three-Dimensional  Space  Marching 


The  consistent  (CSIP)  algorithm  has  also  been  used  for  the  computation 


of  three-dimensional  flow  configurations  with  the  CV/RNS  formulations.  For 


the  problems  that  have  been  considered,  the  axial  flow  direction  (in  an 


appropriate  body-fitted  coordinate  system)  is  discretized  in  a  time-like  or 


marching  character,  and  the  cross-plane  difference  system  is  prescribed  with 


the  familiar  five-point  discretization.  The  global  pressure  relaxation 


algorithm  described  in  previous  discussions  is  directly  applicable  to  these 


problems.  The  flow  over  a  three  dimensional  rectangular  wing  of  zero 


thickness  and  the  flow  past  a  centerbody  of  rectangular  or  elliptical  cross¬ 


section  have  been  considered  to  date.  The  rectangular  elliptical  cross 
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section  ha?  been  examined  with  the  PV/RNS  formulation  and  the  effect  of 
finite  span  on  the  trailing  edge  interaction  has  been  assessed  in  an  earlier 
report.  Current  research  relates  to  variable  camber  or  thickness  effects, 
see  figures  6a, b.  The  CV/RNS  formulation  has  now  been  applied  for  an 
rectangular/elliptical  cylinder  configuration.  The  skin  friction  along  the 
top  and  the  side  symmetry  plane  are  presented  in  figures  6c, d.  The 
composite  velocity/ vorticity  formulation  has  been  modified  to  remove  the 
singular  behavior  that  results  for  zero  incidence  and  zero  cross  flow 
conditions.  The  improved  code  has  been  tested  for  axisymmetric  geometries 
and  is  currently  being  applied  to  the  CV/RNS  computation  of  flow  past  the 
afterbody  configuration  discussed  previously.  These  results  will  appear  in 
a  Ph.D.  dissertation  of  R.  Gordnier  to  be  completed  later  this  year.  A 
detailed  paper  should  be  available  by  early  1989- 

1 "6  Application  of  the  Sparse  Matrix  Direct  Solver  (PS) 

During  the  last  contract  period  an  innovative  approach  to  the  solution 
of  compressible  flow  problems,  based  on  a  direct  sparse  matrix  solver,  had 
been  initiated.  The  solution  algorithm  is  highly  flexible  and  robust.  The 
most  significant  outcome  of  the  research  has  been  the  development  of  a 
domain-decomposition-iterati ve  procedure  based  on  the  DS.  The  computational 
domain  is  decomposed  into  several  overlapping  subdomains.  The  DS  is  then 
used  to  solve  the  problem  in  each  sub-^domain.  Although  a  related  procedure 
has  recently  been  presented  by  others  for  simple  and  highly  elliptic 
equations,  the  present  decomposition  is  found  to  be  more  suitable  for  the 
solution  of  the  Navier-Stokes  equation?  at  high  Reynolds  number.  This 
procedure  honors  the  asymptotic  character  of  the  flow. 


FIG.  6a.  3-D  VARIABLE  CAMBERED  WING 
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The  technique  has  been  applied  to  the  computation  of  flow  past  the 
sine-wave  airfoil  at  high  Reynolds  number  discussed  previously.  Earlier 
calculations  on  coarser  meshes  had  predicted  a  considerable  dependence  of 
the  solution  on  the  grid  distribution;  in  particular,  sensitivity  on  the 
size  and  structure  of  the  separation  bubble  and  possible  laminar  flow 
breakdown  with  NS  and  RNS  solvers.  Recent  calculations  with  the  full  (\jr-oi) 
form  of  the  NavierMStokes  equations  and  with  second  order  central 
differencing  on  a  very  fine  mesh  of  355x101  have  been  possible  using  the 
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domain  decomposition  procedure.  The  skin  friction  at  Re  =  5x10  is  shown  in 
figure  7a.  Solutions  using  seven  domains  converged  in  35~30  iterations. 
Results  for  a  variety  of  Re  are  given  in  figure  7b.  A  paper  on  this 
investigation  will  be  presented  at  the  National  Fluid  Dynamics  Conference  in 
Cincinnati  in  June  1988. 

Additional  issues  relating  to  non-linear  iteration  and  sensitivity  to 
initial  conditions  have  also  been  investigated.  A  paper  on  this  study  will 
be  presented  at  the  International  Conference  on  Computational  Engineering 
Science  to  be  held  in  Atlanta  in  April  1988.  Furthermore,  the  application 
of  the  DS  to  the  computation  of  transonic  inviscid  flow  has  also  been 
completed.  A  modified  Newton  method  and  several  other  techniques  have  been 
investigated  to  explore  the  sensitivity  to  initial  conditions.  Recent 
computations  of  transonic  flows  by  M.  Hafez  using  a  direct  solver  required 
the  application  of  a  continuation  technique  for  the  solution  of  flows  with 
strong  sharp  imbedded  shock.  For  such  cases,  it  has  been  determined  by  the 
present  investigators  that  the  choice  of  initial  conditions  is  critical.  In 
a  majority  of  the  transonic  flow  cases,  the  direct  application  of  Newton's 
method  will  lead  to  divergence.  A  number  of  modifications  have  been  carried 
out  in  order  to  eliminate  this  sensitivity  to  initial  conditions.  The 


solution  for  the  flow  over  a  biconvex  airfoil  at  M  =  0.92  has  been  obtained 
very  efficiently.  This  flow  exhibits  a  complex  fishtail  shock  that  has  been 
computed  with  arbitrary  initial  conditions,  see  figures  7c, d.  These  new 
features  include  the  modification  of  the  Newton's  method  and/or  the  addition 
of  a  transient  artificial  viscosity  that  vanishes  in  the  converged  state. 
More  complete  results  will  appear  in  the  Ph.D.  thesis  of  E.  Bender  later 
this  year.  Several  technical  papers  are  in  progress  and  should  be  completed 
by  early  1989. 

Application  of  the  sparse  matrix  DS  for  three  dimensional  space 
marching  is  also  under  investigation.  The  DS  is  used,  at  each  marcing  step, 
to  solve  for  the  discrete  algebraic  system  in  the  secondary  flow  plane. 
Preliminary  solutions  have  been  obtained  for  flow  past  a  finite  plate  and 
along  an  axial  (90°)  corner.  A  paper  on  the  three-dimensional  wing-bump 
solutions  presented  previously  has  been  presented  at  a  CFD  meeting  in 
Sydney,  Australia  in  August  1987.  A  second  paper  will  appear  in  the  near 
f  ut  u re . 

1 . 7  RNS  Relaxation  and  Flux  Vector  Splitting 

The  relation  between  global  pressure  relaxation  and  a  new  form  of  flux 
vector  splitting  has  recently  been  clarified  by  the  present  investigators. 

A  paper  on  this  subject  is  scheduled  to  appear  in  Computers  and  Fluids.  A 
preliminary  copy  is  appended  herein. 
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1 . 8  Highlights  of  Research  Progress 


1.  An  effective  computational  methodology  based  on  PV  and  CV/RNS  equations 
for  the  solution  of  inviscid  and  viscous  flows  with  strong  pressure 
interaction  has  been  established.  The  procedures  require  reduced 
storage  for  steady  flows  and  are  computationally  efficient  for  steady 
and  transient  flows.  Several  steady  state  and  transient  flow 
applications  have  been  completed  for  airfoil  and  transonic  flows.  An 
alaternating  direction  (predictor-corrector )  algorithm  has  been  applied 
successfully  with  larger  values  of  At  than  allowed  by  NS  solvers  for 
airfoil  interaction  with  shedding. 

2.  A  consistent  version  of  the  coupled  strongly  implicit  procedure  (CSIP) 
for  unsteady  and  three-dimensional  space  marching  flows  has  been 
developed.  A  number  of  PV  and  CV/RNS  solutions  have  been  presented 
with  this  algorithm.  Three-’dimensional  applications  using  the  CV/RNS 
development  have  been  demonstrated  for  an  afterbody  configuration  and 
with  the  PV/RNS  development  for  a  wing  of  variable  camber. 

3.  A  hysteresis  phenomena  for  unsteady  flow  over  airfoils  has  been 
predicted  and  has  now  been  confirmed  with  more  detailed  and  refined 
computations . 

4.  For  separated  flow  computations,  the  strong  influence  of  grid  on  the 
resolution  of  flow  behavior  has  been  established.  A  laminar  flow 
'instability'  has  been  observed  for  fine  meshes  and  large  Reynolds 
number.  Further  RNS  and  NS  analyses  have  demonstrated  that  the 
appearance  of  this  "shock  like"  reattachment  is  probably  physical  and 
not  numerical.  Grid  refinement,  temporal  accuracy  and  large  Re 
asymptotic  theories  have  reinforced  this  result.  Artificial  viscosity 
in  conventional  NS  solvers  suppresses  this  behavior  as  well  as  the  true 
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diffusional  influences  on  the  flow.  This  has  been  shown  independently 
by  F.  Marconi  of  Grumman  and  by  researchers  at  NASA  Langley  for  flow 
over  conical  and  airfoil  geometries. 

5.  A  uni-directional  multigrid  pressure  relaxation  algorithm  has  been 
developed  for  highly  stretched  grids  associated  with  thin  viscous 
layers  and  separated  flows.  Conventional  procedures  fail  in  these 
cases.  This  procedure  is  now  being  considered  for  three-dimensional 
computations . 

6.  The  RNS  philosophy  for  high  supersonic  flows  has  been  confirmed.  Flows 
with  shock-shock  interaction  and  cavities  have  been  successfully 
computed.  A  flux  vector  split  analysis  has  demonstrated  the 
relationship  of  pressure  splitting  to  the  flux  eigenvalues.  This 
procedure  appears  quite  promising  for  hypersonic  flows. 

7.  The  application  of  a  sparse  matrix  direct  solver  (DS)  for  the  solution 
of  separated  flow  and  transonic  flow  problems  has  been  investigated. 

The  application  of  the  DS  can  markedly  improve  the  computational 
efficiency  of  the  solution  procedure.  Overlapping  grid  DS  procedures 
appear  to  be  quite  promising.  A  variety  of  solutions  have  now  been 
obtained. 

8.  Applications  to  three-dimensions  have  been  initiated.  The  influence  of 
moderate  bump  on  a  wing  geometry  has  led  to  the  first  3D  separation 
solutions.  Further  results  have  been  obtained  for  juncture  and 
afterbody  geometries. 
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2.  AFOSR  PUBLICATIONS,  PRESENTATIONS,  ACTIVITIES  AND  INTERACTIONS 


1/87  -  1/88 


Publications  and  Proceedings 

Khosla,  P.K.,  Rubin,  S.G.:  Consistent  Strongly  Implicit  Iterative 
Procedures  for  Two  Dimensional  Unsteady  and  Three  Dimensional  Space 
Marching  Flow  Calculations.  Computers  &  Fluids,  J_5,  9,  PP-  361-373, 
October  1987. 

Lai,  H.T.  and  Khosla,  P.K.:  Global  Pressure  Relaxation  Procedure  for 
Steady  State  Euler  Equations.  Computers  &  Fluids,  1_5,  2,  pp.  215-228, 

July  1987. 

Ramakrishnan ,  S.V.  and  Rubin,  S.G.:  Transient  Flow  Past  a  Finite  Flat 
Plate  at  Incidence.  AIAA  Paper  87-0192,  February,  1987. 

Ramakrishnan,  S.V.  and  Rubin,  S.G.:  Numerical  Solution  of  Unsteady 
Compressible  Reduced  NavieraStokes  Equations.  AIAA  Journal,  2j>,  7,  pp. 
905-91  3,  July  1987. 

Bender,  E.E.  and  Khosla,  P.K.:  Solution  of  Two-Dimensional  Navi er-Stokes 
Equations  Using  Sparse  Matrix  Solvers.  AIAA  Paper  87-603,  February,  1987. 
Himansu,  A.  and  Rubin,  S.G.:  Multigrid  Acceleration  of  a  Relaxation 
Procedure  for  the  RNS  Equations.  AIAA  Paper  87*1195  CP,  June  1987. 
Gordnier,  R.E.  and  Rubin,  S.G.:  Transonic  Flow  Solution  Using  a  Composite 
Velocity  Procedure  for  Potential,  Euler  and  RNS  Equations.  To  appear  in  a 
special  issue  of  Computers  and  Fluids,  July  1988. 

Rubin,  S.G.:  PNS  Methodology  Hypersonic  to  Low  Speed  Flows.  AIAA  Special 
Publication,  June  1987. 

Reddy,  D.R.  and  Rubin,  S.G.:  Consistent  Boundary  Conditions  for  RNS 
Schemes  Applied  to  3D  Internal  Flows.  AIAA  Paper  88-0719,  January  1988. 


V 


10.  Cohen,  R.  and  Khosla,  P.K.:  Three-Dimensional  RNS  Solutions  for  Subsonic 
Unseparated  Flows  Using  Global  Pressure  Relaxation.  Proceedings  2nd 
International  Conference  on  Numerical  Methods  in  Fluid  Mechanics,  Sydney, 
Australia,  August,  1987. 

To  Appear 

11.  Bender,  E.  and  Khosla,  P.K.:  Application  of  Sparse  Matrix  Solvers  and 
Newton's  Method  to  Fluid  Flow  Problems.  Proceedings  of  National  Congress 
on  Fluid  Dynamics,  Cincinnati,  Ohio,  June,  1988. 

12.  Bender,  E.E.  and  Khosla,  P.K.:  A  Modified  Newton’s  method  for  the 
Computation  of  Fluid  Flows.  Proceedings  of  International  Conference  on 
Computational  Engineering  Science,  Atlanta,  Georgia,  April  10-14,  1988. 

13.  Khosla,  P.K.,  Rubin,  S.G.  and  Himansu,  A.:  RNS  Solutions  for  Two 
Dimensional  Transient  and  3-D  Steady  Flow.  Proceedings  of  11th  ICMFD 
Conference,  Will iamsbur g ,  Virginia,  June,  1988. 

14.  Rubin,  S.G.:  RNS  Pressure  Relaxation  and  Flux  Vector  Splitting.  Accepted 
for  Computers  and  Fluids,  1988. 

15.  Lai,  H.T.  and  Khosla,  P.K.,  "Global  Pressure  Relaxation  for  Transonic 
Turbulent  Flows,"  accepted  for  Computers  and  Fluids,  1988. 

16.  Ramakrishnan ,  S.V.  and  Rubin,  S.G.:  ;Numerical  Solutions  of  Unsteady 
Compressible  RNS  Equations.  Proceedings  1st  National  Congress  on  Fluid 
Dynamics,  Cincinnati,  Ohio,  June,  1988. 

Presentations 

1.  Himansu,  A.  and  Rubin,  S.G.:  Multigrid  Acceleration  of  a  Relaxation 

Procedure  for  the  RNS  Equations.  8th  AIAA  Computational  Fluid  Dynamics 
Conference,  Hawaii,  June,  1987. 
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Reddy,  D.R.  and  Rubin,  S.G.:  Consistent  Boundary  Conditions  for  RNS 
Schemes  Applied  to  Three-Dimensional  Internal  Flows,  25th  Aerospace 
Sciences  Meeting,  Reno,  Nevada,  January,  1988. 

Cohen,  R.  and  Khosla,  P.K.:  Three-Dimensional  RNS  Solutions.  Second 
International  Conference  on  Numerical  Methods  in  Fluid  Mechanics,  Sydney, 
Australia,  August,  1987 


Seminars  and  Short  Courses 

1.  Rubin,  S.G.,  Ballistics  Research  Laboratory,  December  11,  1987. 

2.  Rubin,  S.G.,  NASA  Lewis  Research  Center,  May  3,  1987. 

3.  Rubin,  S.G.,  Grumman  Aerospace  Corp.,  March  26,  27,  1987. 

4.  Rubin,  S.G.,  AFWAL,  WPAFB ,  March  4,  1987. 

5.  Khosla,  P.K . ,  NASA  Lewis  Research  Center,  September  10,  1987. 


Student  Presentations 

1.  Himansu,  A.:  "  Multigrid  Acceleration  of  Relaxation  Procedure  for  the  RNS 
Equations".  AIAA  Student  Paper  Competition,  Region  III,  Dayton,  Ohio, 
March  1987.  (Won  Second  Prize) 

2.  Pordal ,  H.S.:  "A  Comparison  of  Conjugate  Gradients  and  Minimum  Residual 
Algorithms  for  the  Acceleration  of  Relaxation  for  Three-Dimensional 
Inviscid  Flow  Past  an  Afterbody  Configuration".  AIAA  Student  Paper 
Competition,  Region  III,  Dayton,  Ohio,  March  1987. 

3.  Himansu,  A.:  "Multi-Grid  Acceleration  of  a  Global  Pressure  Relaxation 
Procedure".  M.S.  Thesis,  January,  1987  (Now  Ph.D.  Student). 

4.  Ramakr i shnan ,  S.V.:  "Numerical  Solution  of  Unsteady  Compressible  Reduced 
Navier-Stokes  Equations".  Ph.D.  Dissertation,  February  1988  (now  at 
Rockwell  International,  Sunnyvale,  California). 
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5. 


Pordal,  H.S.:  "Conjugate  Gradient  Type  Methods  for  Three-Dimensional 


Potential  Flow  Calculations."  M.S.  Thesis,  October  1987  (now  Ph.D. 
Student) . 

Bender,  E. :  "Use  of  Direct  Sparse  Matrix  Solvers  and  Newton's  Iteration 
for  the  Numerical  Solution  of  Fluid  Flows".  Ph.D.  Dissertation,  April, 
1988  (scheduled). 

Lai,  H.T.:  (see  publications),  Post-Doctoral  Fellow  and  Ph.D.  graduate 
(now  at  Sverdrup,  Inc.,  Cleveland,  Ohio). 


B.  Committees  and  Assignments 
Rubin,  S.G. 

AIAA  Session  Chairman  and  Organizing  Committee,  Fluid  &  Plasma  Conference, 
Honolulu,  1987 

AIAA  Session  Chairman,  CFD  Conference,  Honolulu,  1987. 

NASA  Aerospace  Research  &  Technology  Subcommittee  (ARTS),  1986-88. 

NASA  Aerospace  Advisory  Committee,  CFD  Validation  Subcommittee,  1986-87. 
Case  Institute/NASA  Lewis  Institute  for  Computational  Mechanics  in 
Propulsion  (ICOMP)  Advi sor y  Commi ttee  ,  1987-89. 

WPAFB  -  IPA  -  Visiting  Scientist,  1987. 

Editor-in-Chief  -  International  Journal,  Computers  and  Fluids. 

NASA  Review  Committee  of  Aeronautics  Advisory  Committee  -  Supersonic 


Aircraft  Drag  Reduction,  1988-89. 


Khosla,  P.K. 


Editorial  Advisory  Board,  Computers  and  Fluids 
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C .  Interactions 

During  the  period  of  the  current  AFOSR  contract  the  principal 
investigators  have  interacted  technically  with  several  outside  researchers 
and  organizations.  The  CSip  developed  by  the  Pi's  has  been  applied  to 
internal  flow  problems  by  D.  Reddy  of  Sverdrup,  Inc.  in  conjunction  with  the 
NASA  Lewis  Research  Center.  A  paper  on  this  subject  was  presented  at  the 
AIAA  meeting  in  Reno  in  January  1988.  This  procedure  is  also  being 
considered  for  internal  flows  by  R.  Pletcher  of  Iowa  State  University.  The 
RNS  model  and  CSIP  algorithm  are  continuing  to  be  applied  to  hydrodynamics 
problems  by  Raven  and  Hoekstra  at  MARIN,  The  Netherlands  and  to  aerodynamics 
problems  by  D.  Reddy  at  Sverdrup,  M.  Rosenfeld  of  The  Technion  (Israel)  and 
the  NASA  Ames  Research  Center,  M.  Barnett  at  United  Technologies  Research 
Center  (see  AIAA  paper,  Reno  1988)  and  C.  Fletcher  of  the  University  of 
Sydney  in  Australia.  Roger  Cohen,  a  Fulbright  Research  Fellow  from  the 
University  of  Sydney  has  completed  his  work  on  these  techniques  during  a 
two-year  visit  at  the  University  of  Cincinnati.  He  has  developed  the  three 
dimensional  code  for  the  rectangular  wing  with  spanwise  ripples.  The 
results  of  these  studies  were  presented  at  the  International  CFD  Conference 
in  Sydney  during  August  1987.  His  Ph.D.  dissertation  should  be  completed 
during  1988  and  a  second  paper  with  P.  Khosla  is  in  progress. 

During  1987,  Professor  Rubin  spent  a  sabbatical  lave  on  a  joint 
AFOSR( IP A) -Uni versi ty  Resident  Research  Program  at  AFWAL ,  WPAFB.  He  worked 
with  Dr.  J.  Shang  and  Dr.  D.  Rizzetta  of  the  Computational  Aerodynamics 
Group  to  develop  the  RNS  procedure  for  high  speed  three-dimensional  flows. 
The  earlier  AFOSR  supported  work  on  PNS  and  RNS  computations  have  led  to  two 
result?  that  are  now  default  options  in  the  latest  version  of  the  AFWAL  PNS 
code  as  updated  by  Lockheed.  The  minimum  step  size  restrictions 


=  2yM_1/iT  and  the  sublayer  apoproximations  of  Rubin  and  Lin  are  new 

features  of  the  code.  Moreover,  the  boundary  region  or  cross  flow  diffusion 

terms  have  been  added  to  account  for  large  curvature  or  separation  effects. 

The  result  of  this  activity  is  several  papers  on  high-speed  RNS  techniques 
and  flux  vector  splitting  that  will  appear  in  the  near  future. 

The  work  on  the  application  of  sparse-matri x  direct-solvers  has 
generated  significant  interest  around  the  country.  L.  Wigton  at  Boeing, 

D.  Venkata  at  NASA  Ames,  M.  Salas  at  NASA  Langley,  C.  Merkel  at  Penn  State, 
M.  Hafez  at  University  of  California  (Davis)  are  just  a  few  of  the 

researchers  who  are  actively  engaged  in  pursuing  similar  ideas  and  who  have 

requested  further  information  on  our  investigations.  Papers  on  this 
subject,  by  Professor  Khosla  and  E.  Bender,  are  to  be  presented  at 
conferences  in  Atlanta  this  spring  and  in  Cincinnati  this  summer. 

Professor  Rubin  presented  a  two^day  seminar  at  Grumman  Aerospace  on  RNS 
methods.  A  two-’dimensional/axisymmetric  code  has  been  provided  to  the 
Grumman  group  for  evaluation.  Further  discussions  are  ongoing  with 
Dr.  R.  Melnik  relating  to  a  possible  research  project  on  three  dimensional 
aerodynamic  flows.  As  a  result  of  these  discussions,  a  project  on  the 
effects  of  artifical  diffusion  and  cross  flow  (boundary  region)  diffusion 
was  initiated.  This  has  led  to  some  important  conclusions  that  were 
presented  by  Dr.  F.  Marconi  of  Grumman  at  the  1988  Reno  meeting. 

Professor  A.  Polak  is  continuing  the  application  of  the  RNS  code  for 
the  investigation  of  roughness  effects  at  supersonic  speeds.  The  geometry 
under  consideration  is  a  hollow  cylinder  at  M  =3  for  aircraft 

oo 

considerations.  Professor  P.  Disimile  is  investigating  the  roughness 
problem  in  collaboration  with  the  experimental  aerodynamics  group  at  WPAFB. 
He  has  a  mini-grant  to  aid  in  this  project.  The  RNS  code  as  modified  by 
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Professor  Polak,  will  provide  the  computational  support  for  this  work.  In 
addition,  W.  Sturek  at  BRL  has  expressed  particular  interest  for  application 
to  roughness  problems  in  ballistics.  A  seminar  by  S.G.  Rubin  was  presented 
at  BRL  this  winter.  Additional  computations  by  several  Ph.D.  students  have 
demonstrated  the  strong  convergence  properties  of  the  RNS  formulation  for 
high  supersonic  Mach  numbers  and  strong  pressure  interactions,  including 
shock  waves  and  separation. 

A  new  high  Mach  number  RNS  initiative  with  members  of  the  staff  at 
Sverdrup,  Inc.  is  under  discussion.  Professor  Rubin  will  visit  ICOMP  and 
NASA  Lewis  and  present  a  seminar  at  Sverdrup  this  spring.  Other 
interactions  with  AFWAL,  General  Dynamaics  and  Grumman,  as  they  relate  to 
hypersonic  aerodynamics,  are  under  discussion. 


-.N3  PRESSURE  RELAXATION  AND  FLUX  VECTOR  SPLITTING 

S.G.  Rubin 

Department  of  Aerospace  Engineering  i  Engineering  Me cranio? 
University  of  Cincinnati 
Cincinnati,  OK  ^5221 


1 .  I rtroduct ion 

The  present  author  and  co-workers  have  previously  presented  composite 

primitive  variable'1  ^  and  velocity/ potential '  ^  formulations  for  the 
computation  of  viscous/ invisci 1* interacting  flows.  These  procedures  have 
been  designed  to  emulate  large  Reynolds  number  (Re)  asymptotic  behavior, 
but,  at  the  same  time,  to  allow  for  the  efficient  solution  of  a  single 
composite  system  of  equations.  This  methodology  is  intermediate  to  that  of 
full  time-dependent  Navier-Stokes  solvers  and  to  that  of  interacting  or 
inverse  matcr.ed  boundary  layer- inviscid  solvers.  The  resulting  system  of 
reduced  Navier-Stokes  or  RNS  equations  is  a  composite  of  the  full  Euler  and 
second-order  boundary  layer  systems.  The  neglected  Navier-Stokes  or 
diffusive  terms  are  higher-order  in  Re  for  appropriate  ’streamline’ 
coordinates.  Although,  these  terms  can  be  retained  through  an  explicit 
def erred-corrector ,  it  must  be  emphasized  that  the  boundary  conditions  and 
the  discrete  approximation  to  the  differential  equations  are  dictated  solely 


by  the  form  of  the  lowest-order  implicit  RNS  operator. 

The  RNS  model,  which  is  V3iid  throughout  the  entire  Mach  number  (M) 

range,  represents  an  enhancement  of  the  hypersonic  merged  layer^’10^  and 

thin  layer  approximations^  ,  wherein  the  axial  (x)  or  streamwise  pressure 

gradient  (p  )  is  prescribed,  or  of  the  PNS  model  1  2’ 1  ^  ,  where  p  is 


approximated  with  the  truncated  form  wp^.  The  parameter  oo(M)  is  a  function 


of  the  local  Mach  number  and  is  such  that  0  <  ui  <  .j  <  '  ,  wr.ere  1  :  =  2, 

a)  (1)  =  1.  These  system?  of  equations  are  all  '  hyperbcl  1  o'  ,  i.e.,  exhibit 

initia.- value  character  in  the  axial  or  x  direction.  Forward  marching  or 
backward  (upwind)  differencing  is  applied  in  order  to  obtain  a  solution  in  a 
single  x  sweep.  These  methods  are  inadequate  however,  when  axial  flow 
separation,  transonic  regions  with  strong  shock  waves,  pressure- viscous 
layer  interaction,  large  axial  curvature  or  subsonic  to  moderate  supersonic 
free  stream?  are  present. 

With  appropriate  treatment  of  the  elliptic  character  of  the  equations, 
the  RMS  approximation,  which  allows  for  the  full  implicit  pressure  gradient 
(p  ),  has  been  shown  to  accurately  represent  axial  reversed  flow  and  strong 

pressure- viscous  interaction  and  to  provide  a  more  efficient  method  than 
full  time-dependent  Navier-Stokes  solvers  for  this  important  class  of  flow 
problems.  Although  the  upwind  character  of  the  solution  procedure  is 
retained,  the  RNS  solver  is  no  longer  of  initial  value  type  and  therefore 
the  unwinding,  or  more  appropri ately ,  flux-vector  splitting,  must  reflect 
the  pressure  (or  potential)  and/or  reversed  flow  boundary  value  character. 

A  global  pressure  relaxation  or  multi-sweep  procedure  for  the  RNS 

( 1-4) 

model  has  been  presented  .  Viscous  RNS  and  inviscid  (Re  *  «)  or  Euler 

solutions  have  been  obtained  for  a  variety  of  flows  with  strong  pressure 

(10)  ( n  i u )  (23) 

interaction  ’  ~  ,  with  captured  shocks  ’  and  with  flow  reversal  ’ 

( 1  -  4 )  (14) 

Both  steady  and  transient  algorithm?  have  been  developed  and  a  uni¬ 

directional  (x)  multi-grid  accel erator ^ 1 v ^  has  been  effective  for  improving 
the  convergence  rate  of  the  relaxation  procedure. 


% 

J£ 


trow 


The  RNS  discretization  reflects  first  or  second-order  upwindir.g  only  ir. 


the  axial  or  x  direction.  The  normal  (y)  gradients  are  approximated  with 
'central'  discreti zations .  These  are  two-point  or  trapezoidal  for  first- 


order  terms,  e.g.  v 


p^  in  the  continuity  and  normal  momentum  equations, 


respectively.  This  allows  for  consistent  non-ref leeti v?  boundary  conditions 


and  reduces  numerically  induced  oscillations  that  may  appear  with  full 
three-point  central  discretizations  for  these  terms.  This  also  allows  for 
more  fficient  application  of  far-field  boundary  conditions  and  more 
accurately  represents  the  behavior  of  thin  viscous  layers,  see  references 
1 4)  for  further  details. 

Cther  investigators  concerned  with  Euler,  thin  layer  Navier-Stokes  or 
full  Navier-Stokes  equations  have  also  presented  algorithms  that  are 
designed  to  relieve  the  'artificial-viscosity'  difficulties  associated  with 
centr al-dif f erence  discretizations.  These  are  'upwind'  approximations  that 
are  associated  with  the  movement  of  the  physical  forward  and  backward  moving 
waves,  i.e.,  flux-difference  splitting,  or  with  the  movement  of  the  discrete 
forward  and  backward  moving  particles,  i.e.,  flux-vector  splitting. 

In  previous  PNS  or  RN3  investigations,  the  choice  of 
pressure/convective  upwinding  has  resulted  from  characteristi c  or  stability 
considerations.  However,  these  results  should  also  be  derivable  from  a  form 
of  flux-vector  splitting,  since  the  characteristi c  analysis  is  simply  a 
reflection  of  the  eigenvalue  or  flux  behavior.  Therefore,  in  the  present 
note,  the  pressure-velocity  decomposition,  that  has  been  applied  for  the 
initial  value  PNS  or  boundary  value  RNS  methodology,  is  re-examined  from  the 


point  of  view  of  flux  vector  splitting.  The  results  represent  a  variation 


For  the  purpose  of  the  present  analysis  we  consider  t.ne  one- di sens lona 
Euler  system  given  by 

G_  «■  [ E ( G ) }  =  G„  *  (E+(Q)}  +  {e'(G)}v  =  0  ,  (la) 

W  X  Ur  A  X 

T 

where  G  is  the  solution  vector  Q  =  [p,  ou,  pe]  ,  e  is  the  specific  interna 

+■  — 

energy,  p  the  density  and  u  the  velocty.  E(Q)  =  E  (Q)  +  E  (G)  and  is  such 

that  all  eigenvalues  A"  of  dE+/dQ,  dE  /dG  are  non-negative  and  non-posi ti ve . 
respecti vely . 

In  order  to  reflect  the  initial  value  character isti c  behavior  for 

supersonic  flows,  we  also  require  that  E  (Q)  =  E(Q),  E  (Q)  =  0  for 

2 

M  =  u/a  >  1,  where  a  is  the  sound  speed  (a  =  Yp/p  for  a  perfect  gas),  Y  is 

the  ratio  of  specific  heats  and  p  is  the  pressure  (p  =  (Y-I)pe).  Other 

desirable  propert  i  es  v 1  ^  ^  are  that  E^Q),  dE±(Q)/dQ  and  A1  be  continuous  for 

ail  M,  or  for  all  (— - )  for  incompressible  flow,  and  that  dEVdQ  have  at 

ref 

least  one  zero  eigenvalue  for  | M  j  <  1.  This  latter  condition  allows  for  a 

0  7) 

steady  state  shock  structure  with  a  maximum  of  two  interior  zones  and 
therefore  results  in  more  accurate  shock  capturing.  In  the  present 

(17) 

analysis,  we  have  abandoned  both  the  uniqueness  condition  that  E(Q)  is 

(17) 

the  lowest  degree  polynomial  in  M  and  the  symmetry  condition  for  M  ■*  -M 
The  objective  herein  is  to  identify  the  flux-splitting  as  strongly  as 
possible  with  the  respective  convective  (axial)  and  acoustic  disturbances  in 
order  to  maximize  the  ’marching'  or  relaxation  bias  of  the  RNS  formulation. 
This  is  in  keeping  with  the  methodology  of  interacting  boundary  layer 
calculations,  in  which  the  axial  pressure  gradient  provides  the  only 


I 

I 
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boundary  value  contribution  for  attached  flows.  This  is  also  true  for  the 
Euler  system,  where  the  acoustic  upstream  influence  is  also  pressure  based. 
For  regions  of  reversed  flow,  the  convective  fluxes  must  be  appropriately 
upwir.ced  in  orde”  to  reflect  the  direction  of  the  particle  flow.  The 

(17) 

specification  of  the  acoustic  and  mass  flux-splitting  is  non-unique 

Two  new  flux  vector  splittings  that  satisfy  many  of  the  properties  on 

E~  delineated  previously  are  introduced  here.  In  one  model,  the  continuity 

condition  on  E-  at  separation  and  reattachment  points  is  no  longer 
satisfied.  For  a  second  more  satisfactory  model,  ail  of  the  properties 
including  continuity  are  satisfied;  however,  the  additional  requirement  that 

E  (Q)  =  E(5),  E+(Q)  =  0  for  M  £  -1  is  not  enforced.  Since  we  are  not 
concerned  here  with  reversed  supersonic  regions,  this  additional  condition 
does  not  appear  to  be  a  severe  limitation.  Morever,  it  is  possible  to 
'blend'  the  two  forms  of  flux  splitting  in  order  to  move  smoothly  from  one 
to  the  other  and  thereb:  satisfy  all  conditions.  This  would  lead  to  a 

'central'  difference  approximation  over  a  limited  range  where  -  |  M  j  £  M  £  0 

or  -  |  u  |  £  u  £  0 . 

In  order  to  evaluate  the  pressure/convecti ve  form  of  flux-vector 
splitting,  we  consider  the  simplified,  partially  conservative,  Euler 
system  (1).  The  differential  system  is  given  by 


♦  (dE  /dQ)Qx  +  (dE  /d£ )  Qx  =  0 


or  in  scalar  form 


pt  *  (pu)x  =  0 


(pu)  -  u ‘p  ♦  2u(pu)  ♦  w(7-1)(pe)  +  ( 1  -u) )  ( Y- 1  )  (  pe ) 

t  X  X  X  X 


( pe), 


u(pe)  *  e(pu)  -  euo  *  ^  (pu)  -  p  =  0 
X  X  X  p  xpx 


(ie) 


where  ui  =  w'M). 

In  the  HNS  methodology  for  attached  flows,  all  spatial  gradients,  save 
the  pressure  gradient  (*-u>)(Y-l  )  (pe)  term  in  (Id),  are  associated  with 

E  (G).  Therefore 


0 

1 

0 

1 

0 

G 

0 

2 

-u 

2u 

o)(Y*1  ) 

dE~ 

’  dQ 

0 

0 

(1-w)(Y-l ) 

-uh 

h 

u 

0 

0 

0 

(2) 

where  h  =  e  *  p/p  and  a  =  (Y-1)h  for  a  perfect  gas.  The  values  of  E  i3), 
E  (Q)  can  be  approximated  by  the  following  generalized  expressions: 


E * ( Q )  =  [pu,  puu 


w  ( x )  _  u  ( x ) 

(Y-1)t>pe  -  /  pe  d£],  pue  +  J 
0 


puvx . 


r  pd E,  +  I 
0 


w(x)  .  u(x) 

E  (G)  =  [o,  puu  +  ( Y-1 )[ ( 1 -w)pe  +  /  pe  d£],  oue  +  /  pd- 

0  0 


h  dtj 

(3a) 

pu(x) 

/  hdE]1 

0 

(3b) 


where  u  =  (u,0)  and  u  =  -(-u,0)  .  Note  that  the  model  equations  are 

’  max  max 

not  in  full  conservation  form  and  this  leads  to  the  integral  flux 

contribution  in  E±(Q).  Also  note  that  the  contribution  in  E"(Q)  from  the 

2 

pressure  gradient  p^  or  (Y-1)(pe)x  in  (Id)  is  not  given  by  pu  +  ;jp  and  (1- 
w)p,  respectively.  The  parameter  w  is  functionally  dependent  on  M  and  has  a 
discontinuity  in  slope  at  M=l.  However,  with  the  definitions  (3)  for 
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a)  contn  but-ions  will  not  appear  in  dn-(Q  )/ci .  These  term?  are  not  present 
in  dE(Q)/dQ.  Therefore,  for  the  RNS  splitting  given  by  (2,  3'.  the 
continuity  of  dE~(Q)/dQ  is  maintained  at  sonic  points  to  insure  smooth  shocK 

4- 

capturing.  The  present  form  of  E~(3)  leads  directly  to  the  RNS  pressure 
decomposition  that  has  been  applied  in  earlier  analysis.  The 

discretizations  for  (E~(Q) can  now  be  given  by 


[Vi  •  -  1  P‘»C/<vxi-i) 


f cr  E  (Q),  and 


:(1'^i  +  l  )pi  +  1  "  O-oii  )Pj,  -  /  P  <U]/(xi  +  rx.) 


for  E  (Q) 


For  J  p  d£  -  (u>.  -  u).  ,  )  p  and  p  e  [p.  .  ,d.],  we  obtain  for 

l  i-1  *i- 1  •  l 


(E  <3>}x. 


0)  .  +  CJ  -  , 

1  1-1 


(p.-p. 

1  1- 


for  p  = 


i  'i-1 


and  for  ( E  (3 ) }  , 


1  ->  U)  .  , 

1  ♦  1 


( p .  - p .  )  for  p 

i  +  l  l  K 


(’■  V 


^hese  are  precisely  tne  discrete  approximation?  specified  in  previous  F.NS 
computations . 

The  eigenvalue  (X-)  analysis  for  the  matrices  (2)  leads  to  the 
following  for  u  >  0; 

*  ,  1/2  1/2  .  , 
a,  =  vu  +  aw  ,  u  -  aui  ,  u)  ,  (5a) 

'  t  *-  t  J 

ail  of  which  are  non- negative  for  u  >  0,  and 

X"  -  =  (0,  0,  0)  ,  (5b) 

1  I  ^  f  J 

ail  of  which  are  zero;  therefore,  in  region?  of  attached  flow  upstream 
influence  is  of  purely  diffusive  character.  These  results  require  that  the 


function  m  is  given  by  0  <  cu  <  uj  =  M  .  Fcr 

—  —  m 


X1  2  ,  =  ( 2u ,  0,  u)  (5c) 

and  therefore  d£  (2)/dQ  has  one  zero  eigenvalue  for  all  |M  j  <  1. 

For  separated  flows  (u  <  0),  one  choice  of  flux  splitting  is  to  reverse 


the  roles  of  all  fluxes  so  that  E  +  (Q)  -*•  E+(Q)  and  hence  \t  *■  X+.  As  noted 
previously,  this  satisfies  all  of  the  positivity  (negativity)  conditions. 

Although,  the  eigenvalues  of  3E_(Q)/cQ  are  continuous,  both  E~(Q)  and 

dE“(Q)/dQ)  are  discontinuous  when  u  change?  sign.  Althougn  this  splitting 

(14) 

has  been  applied  successfully  for  several  time  dependent  computations'  , 
it  is  sensitive  to  the  temporal  increment  and  marginally  stable  fcr  Newton 
linearized  steady  state  relaxation.  An  alternate  and  more  desirable  form  is 
ootained  with  the  fluxes  (2,3).  but  with  w  =  0  for  all  u  <_  0.  The  flux 
expression?  (3)  for  u  <  0  become 

pu(x) 

E+(Q)  =  [pu,  0,  /  h  d£]‘  (6a) 

0 


»>  V  W  V  V  V  V/W  v  V  V'W^W 


so  that  for  u  <  0 


o 

o 

0  0  G 

0  G  0 

dE  (C) 

-u“  2u  Y- 1 

’  cQ 

o 

o 

-  uh  0  u 

In  this  f  orm  E  (Q ) 

and 

dE+(S )/dQ 

are 

continuous  at 

u  =  0,  as  well  as  at  M= 

as  are  the  eiger.va 

1  ues 

X~,  which 

f  or 

u  <  0  are  all 

zero  or  non- pos i ti ve , 

res  pecti vel y : 

x;  =  (o, 

0, 

0) 

(7a) 

=  ( n 

A1,2,3  "  ’ 

2u , 

u ) 

(7b) 

Identical  results  for  the  eigenvalues  X—  are  obtained  with  equations 
(1)  given  in  full  conservation  form  in  terms  of  entropy,  or  in  terms  of 

2 

total  energy  (e^  =  e  +  u  /2),  if  the  up^  work  term  in  the  resulting  energy 
equation  is  split  in  the  same  manner  as  is  p^  in  the  momentum  equation.  If 
the  up^  term  in  the  total  energy  equation  is  unsplit,  i.e.,  constant 
stagnation  enthalpy  or  fully  backward  differenced,  an  alternate  form  of  flux 
splitting  results  and  the  eigenvalues  X-  for  u  >  0  become 

X*  -  [u,  5  (2  +  (Y-l)O-ui)  ±  3({CY-1  )M(1-Ui)}2  +  U^1'2]  (  8  a  ) 

'  t  <-  I  J 

and 

x~  =  [o,  ( i-u)u( y-i ) ,  o]  ict ; 


•i 

y 


\ 


V  Wj 


where  again  0  <  to  <  w  ,  but  now  ^ 

e  -  -  m’  m 


1  +  ( 1  )M" 


—  and  x7  =  0  for  M  <  1 


for  ui  =  .  This  is  a  slightly  different  condition  than  obtained  in  (5)  and 

reflects  the  modification  of  the  up^  flux  splitting.  This  condition  was 

f  1 2 ) 

obtained  previously  by  Vigneron  et  al .  ‘  ,  where  it  was  also  assumed  that 

E  (Q)  =  0,  so  that  an  initial  value  or  fully  backward  differenced  or  forward 
marched  PNS  solver  results.  Of  course,  as  noted  previously,  this  truncated 
method  is  in  considerable  error  for  low  speeds,  for  separated  flows  and  for 
attached  flows,  where  upstream  pressure  interaction  is  non-negl i gi ble ,  e.g. 
strong  shocks,  large  curvature,  surface  discontinuities,  etc.  The  RN’S 

formulation  retains  the  full  E  (Q)  contribution  for  all  forms  of  the 
differential  approximation.  The  diffusion  terms  that  are  omitted  from  the 
full  Navier-Stokes  system  can  also  be  recovered  with  a  global  deferred- 

( 1  U ) 

corrector,  when  such  effects  are  not  dominant 

For  reversed  flows,  full  shifting  of  the  E-(Q)  terms  is  also  applicable 
fcr  the  total  energy  system;  however,  once  again  the  alternative  splitting 
wherein  for  u  <  0  we  specify  only  u  =  0,  is  preferable.  This  leads  to  the 

desired  continuity  of  E±(Q)  and  dE±(Q)/dQ.  For  the  total  energy  forms  of 

the  equations  the  eigenvalues  X±  for  u  <  0  are, 


X  ,  =  ( 0 ,  -  (  Y- 1  )  u  ,  0 ) 

*  t  »  -• 


X,  0  ,  =  (0,  Yu,  2u: 


When  the  total  energy  is  one  of  the  dependent  variables,  an  additional 
condition  is  required  to  maintain  the  full  upwind  form.  Since 
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y-  ]  ? 

P  =  (Y-l)pe  -  —r—  pu“,  the  (pu)  flux  arising  from  p  also  change?  sign 

I  ^  X  X 

with  u;  therefore,  these  terms  which  appear  in  the  momentum  and  energy 

equations  must  also  be  given  in  u,  u  form.  The  results  (9)  are  then 
recovered.  This  additional  condition  does  r.ct  arise  when  the  temperature  or 
the  specific  internal  energy  is  one  of  the  dependent  variables;  the  results 
(7)  are  then  recovered. 

3 .  Stability 

(19) 

It  can  be  shown  that  both  the  full  shift  and  u  =  0  shift  flux- 
splittings  are  unconditionally  stable  for  steady  state  iine  relaxation. 

These  methods  are  sensitive,  however,  to  the  assumed  initial  conditions  for 
pressure,  etc.  that  are  required  for  steady-state  calculations ,  for  terms 
-  (19) 

entering  from  E  (Q).  It  has  been  shown  that  this  sensitivity  can  be 
controlled  by  (i)  an  improved  initial  guess,  e.g.,  several  inviscid 
relaxation  steps,  (ii)  local  non-linear  iteration  at  each  axial  location  or 
(iii)  underrelaxation  of  the  velocities  in  separated  regions,  by  the 
inclusion  of  temporal  damping,  i.e.,  finite  At,  in  the  early  stages  of  the 
relaxation  process.  After  a  few  steps,  e.g.,  5  to  20,  full  relaxation 
without  local  iteration  can  be  restored.  It  can  be  shown  that  the 

discreti zation  for  the  elliptic  pressure  term  arising  from  E  (I)  correspond? 
to  complete  overrelaxation .  While  tnis  is  desirable  close  to  convergence, 
it  is  undesirable  in  the  initial  stage?  of  the  relaxation  process .  This 
difficulty  can  be  tempered  by  a  pressure  correction  term  developed  by 
(18) 

Israeli  .  The  correction,  which  vanishes  in  the  steady  state  introduces 


pressure  underrelaxation  and  also  residual 


nothing  for  multi-grid 


transfer?.  A  combination  of  the  pressure  smoothing  with  full  overrelaxation 
has  provided  a  useful  approach  for  the  semi-coarsenir.g  or  uni-directional 

(15  19) 

multi-grid  accelerator  ’  . 

9.  Conclusion 

Global  pressure  relaxation  for  a  pressure  split  form  cf  the  Euler  or 
reduced  Navier-Stokes  (RNS)  equations  has  been  shown  to  be  equivalent  to  a 
form  of  flux-vector  splitting  that  satisfies  the  major  eigenvalue  and 
continuity  constraints  on  the  fluxes  and  flux  derivatives.  The  upwinding  is 
applied  only  in  the  axial  or  'mass-flow'  direction.  In  keeping  with  the 
asymptotic  form  of  the  RNS  system,  two-point  or  trapezoidal  discreti zation 
is  used  for  appropriate  normal  gradients  in  order  to  allow  for  the  accurate 
evaluation  of  shear  layers  and  a  consistent  specification  of  far  field 
boundary  conditions.  If  the  pressure  gradient  parameter  u>(M)  is  given  by 

its  maximum  value  in  subsonic  regions,  one  eigenvalue  of  X+  is  always  zero 
and  therefore  shock  resolution  is  improved.  For  regions  of  reversed  flow 
convective  upwinding  is  combined  with  the  condition  uj  =  0.  This  insures 
that  the  fluxes,  flux  derivatives  and  eigenvalues  remain  continuous 
throughout  the  flow.  This  form  of  flux  vector  splitting  is  specifically 
designed  to  maintain  a  bias  in  the  direction  of  the  convective  fluxes  and 

(17) 

therefore  abandons  the  (±)  symmetry  of  the  earlier  forms  of  flux  vector 
splitting.  The  upwinding  leads  to  a  relaxation  method  that  is  solely 
acoustic  driven  throughout  sucsonic  regions,  but  also  includes  convective 
relaxation  in  regions  of  reversed  flow. 
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